#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define apd   (8)
#define nat   (apd*apd*apd) //attention le nombre d'atome dans la structure est de nat (0->nat-1)
#define dt   (0.01)
#define dt2  (0.5*dt)   // nat


typedef float Real;

Real alea()
{
	return (Real) ( (0.5 + (double) rand()) / ( 1.0 + (double)RAND_MAX) ); 
}

Real sig     =  0.65;
Real Epsilon =  0.1;

typedef struct 
{
	Real x;        
	Real y;
	Real z;
	Real w;
} point;

point dim = {5,5,5,0};


void pbc( point *r )
{
	const point dim2 = {dim.x/2.0,dim.y/2.0,dim.z/2.0,0};
	
	while(r->x > dim2.x) { r->x -= dim.x; }
	while(r->x <-dim2.x) { r->x += dim.x; }
	
	while(r->y > dim2.y) { r->y -= dim.y; } 
	while(r->y <-dim2.y) { r->y += dim.y; }
	
	while(r->z > dim2.z){ r->z -= dim.z; }
	while(r->z <-dim2.z){ r->z += dim.z; }
}

Real eval_forces( point *force, point *acc, const point *pos, Real *D_ptr )
{
	Real V = 0;
	Real D = 0;
	int i,j;
	
	// set forces to zeros
	memset( &force[0], 0, nat * sizeof(point) );
	
	// accumulate
	for( i=0; i < nat; i++)
	{
		const point *r_i = &pos[i];
		point       *f_i = &force[i];
		//const point dim2 = {dim.x/2.0,dim.y/2.0,dim.z/2.0,0};
		
		for( j=i+1; j < nat; j++)
		{
			Real         V_ij = 0;          //a voir
			const point *r_j = &pos[j];
			point       *f_j = &force[j];
			
			point  r_ij  = { 
				r_j->x - r_i->x,
				r_j->y - r_i->y,
				r_j->z - r_i->z, 0};
			

			pbc( & r_ij );
			
			const Real     r2  = r_ij.x * r_ij.x + r_ij.y * r_ij.y + r_ij.z * r_ij.z;
			const Real	   r4  = r2*r2;
			const Real     r8  = r4*r4;
			const Real    r14  = r8*r4*r2;
			const Real  Sigma6 = sig*sig*sig*sig*sig*sig;
			const Real Sigma12 = Sigma6*Sigma6;
			const Real       A = 48 * Epsilon * Sigma12;
			const Real       B = 24 * Epsilon * Sigma6; 
			
			point         f_ij = { 0 };
			
			const Real   lj  = (sig*sig) / r2; //q^2
			const Real   lj3 = lj*lj*lj;       //q^6
			
			//-- potentiel
			V_ij = 4.0 *Epsilon*( (lj3*lj3) -  lj3 );
			V += V_ij;  //4*eps*((q^6)^2-q^6)
			
			Real Pot = A / r14 - B / r8;
			
			//printf( "r2_%d,%d= %g | V= %g\n", i, j, r2, V);
			
			//-- f_ij
			f_ij.x  = Pot * r_ij.x;
			f_ij.y  = Pot * r_ij.y;
			f_ij.z  = Pot * r_ij.z;
			
			f_i->x -= f_ij.x;
			f_i->y -= f_ij.y;
			f_i->z -= f_ij.z;
			
			f_j->x += f_ij.x;
			f_j->y += f_ij.y;
			f_j->z += f_ij.z;
			
			D += f_ij.x * r_ij.x + f_ij.y * r_ij.y + f_ij.z * r_ij.z;
			
		}
		
		//D +=  r_i->x * f_i->x + r_i->y * f_i->y + r_i->z * f_i->z;
		
	}
	
	*D_ptr = D;
	
	//compute acceleration
	for( i=0; i < nat; i++)
	{
		const point *r = &pos[i];
		point       *a = &acc[i];
		const point *f = &force[i];
		
		const Real   one_over_mass = r->w;
		
		a->x = f->x * one_over_mass;
		a->y = f->y * one_over_mass;
		a->z = f->z * one_over_mass;
		//fprintf( stderr, "a=[%g %g %g]\n", a->x, a->y, a->z);
	}
	
	// return potential
	//fprintf( stderr, "V= %g\n", V);
	return V;	
}

void modif_v_et_r( point *pos, point *vel, const point *acc)
{
	int i;
	for( i=0; i < nat; i++)
	{
		point *r = &pos[i];
		point *v = &vel[i];
		const point *a = &acc[i]; 
		
		v->x += dt2 * a->x;
		v->y += dt2 * a->y;     //new vel
		v->z += dt2 * a->z;
		
		r->x += dt * v->x;
		r->y += dt * v->y;		//new pos
		r->z += dt * v->z;	
		
		pbc( r );               // PBC, optional
		
	}
}

void modif_v_only(point *vel, const point *acc)
{
	int i;
	for(i=0;i<nat;i++)
	{
		point       *v=&vel[i];
		const point *a=&acc[i];
		v->x += dt2 * a->x;
		v->y += dt2 * a->y;
		v->z += dt2 * a->z;
	}
}

Real eval_K( const point *vel )
{
	int i;
	Real K = 0;
	for(i=0;i<nat;i++)
	{
		const point *v  = &vel[i];
		Real         v2 = (v->x * v->x) + (v->y * v->y) + (v->z * v->z);
		
		K += v->w  * v2;
	}
	
	//printf("K= %g\n", 0.5 * K );
	return 0.5 * K;
}


int main(){
	//on déclare les différentes variables et on les initialise
	FILE *fichier = NULL;
	FILE *fichier_temp_press = NULL;
	
	Real  Vol = dim.x*dim.y*dim.z;
	
	int  i,j,k,npas=0,maxpas=5000,label=0,label2=0,kb=1;
	float  l,m,n;
	int  regular = 1;
	
	Real K=0, V=0, T=0, P=0, D=0,T1=0.27,T2=0,K_th = (3.0/2.0)*nat*kb*T1;
	Real mass    = 1.0;
	Real   a=0;
	Real   b=0;
	
	srand( time(NULL) );
	
	//point pos[nat], vit[nat], acc[nat], force[nat];
	point *pos   = (point *)calloc( nat, sizeof(point) ); 
	point *vit   = (point *)calloc( nat, sizeof(point) ); 
	point *acc   = (point *)calloc( nat, sizeof(point) ); 
	point *force = (point *)calloc( nat, sizeof(point) ); 
	
	// initialisation des part	et force
	if( regular )
	{
		Real dL = pow( Vol/nat, 1.0/3.0 );   //pow(a,b) 
		Real ecart = 0.1 * dL;
		printf("dL=%g\n", dL);
		for(i=0;i<apd;i++)
		{
			for(j=0;j<apd;j++)			
			{
				for(k=0;k<apd;k++)
				{
					pos[label].x= (i*dL) + (0.5 - alea() ) * ecart;
					pos[label].y= (j*dL) + (0.5 - alea() ) * ecart;
					pos[label].z= (k*dL) + (0.5 - alea() ) * ecart;
					pos[label].w = 1.0 / mass;
					
					vit[label].x= 0;
					vit[label].y= 0;
					vit[label].z= 0;
					vit[label].w= mass;
					
					acc[label].x=0;
					acc[label].y=0;
					acc[label].z=0;
					
					++label;
				}
			}
		}
	}
	
	
	else
	{
		Real dL = pow( Vol/nat, 1.0/3.0 );   //pow(a,b) 
		//Real ecart = 0.1 * dL;
		for(i=0;i<apd;i++)
		{
			for(j=0;j<apd;j++)			
			{
				for(k=0;k<apd;k++)
				{								
					pos[label].x= (i*dL);
					pos[label].y= (j*dL);
					pos[label].z= (k*dL);
					pos[label].w = 1.0 / mass;
				
					vit[label].x= 0;
					vit[label].y= 0;
					vit[label].z= 0;
					vit[label].w= mass;
					
					acc[label].x=0;
					acc[label].y=0;
					acc[label].z=0;
					
					label++;
				
								
								//else
								//{
								//	label++;
								//	continue;
								//}

				}
			}			
		}	

		for(l=0;l<apd-1;l++)
		{	
			for(m=0.5;m<apd-1;m++)			
			{
				for(n=0.5;n<apd-1;n++)
				{
					pos[label].x= (l*dL);
					pos[label].y= (m*dL);
					pos[label].z= (n*dL);
					pos[label].w = 1.0 / mass;
				
					vit[label].x= 0;
					vit[label].y= 0;
					vit[label].z= 0;
					vit[label].w= mass;
					
					acc[label].x=0;
					acc[label].y=0;
					acc[label].z=0;
					label++;
		//printf("%f %f %f %d\n",n,m,l,label);
				}
			}
		}	
		for(l=0.5;l<apd;l++)
		{	
			for(m=0;m<apd;m++)			
			{
				for(n=0.5;n<apd;n++)
				{
					pos[label].x= (l*dL);
					pos[label].y= (m*dL);
					pos[label].z= (n*dL);
					pos[label].w = 1.0 / mass;
				
					vit[label].x= 0;
					vit[label].y= 0;
					vit[label].z= 0;
					vit[label].w= mass;
					
					acc[label].x=0;
					acc[label].y=0;
					acc[label].z=0;
					label++;	
				}
			}
		}
		for(l=0.5;l<apd;l++)
		{	
			for(m=0.5;m<apd;m++)			
			{
				for(n=0;n<apd;n++)
				{
					pos[label].x= (l*dL);
					pos[label].y= (m*dL);
					pos[label].z= (n*dL);
					pos[label].w = 1.0 / mass;
				
					vit[label].x= 0;
					vit[label].y= 0;
					vit[label].z= 0;
					vit[label].w= mass;
					
					acc[label].x=0;
					acc[label].y=0;
					acc[label].z=0;
					label++;
				}
			}
		}
	}
	
	//printf( "RAND_MAX= %d\n", RAND_MAX );
	Real C = 2.0/(3.0 * nat * kb);

	//boucle principale
	for( i=0; i < nat; ++i )
	{
		pbc( &pos[i] );
	}
	
	// position initiale
	fichier = fopen("atoms.xyz", "wb");
	if(fichier)
	{
		fprintf( fichier, "%d\n", nat );
		fprintf( fichier, "\n" );
		for( i=0; i < nat; ++i )
		{
			const point *r = &pos[i];
			fprintf(fichier,"H %lf %lf %lf\n",r->x,r->y,r->z);
		}
		fclose(fichier);
	}
	
	
	
	
	// eval_forces(...)
	V = eval_forces( &force[0], acc, pos, &D );
	K = eval_K( vit );
	printf("E=%10g  V=%10g K=%10g\n", V+K, V, K);
	
	
	
	fichier = fopen("nrj.txt", "wb" );
	if( fichier )
	{
		fprintf( fichier, "0 %g %g %g\n", V, K, V+K);
		fclose( fichier );
	}
	
	fichier_temp_press = fopen("temppress.txt","wb");
	if(fichier_temp_press)
	{
		fprintf(fichier_temp_press,"%d  %lf  %lf\n",npas,T,P);   
		fclose(fichier_temp_press);
	}

	for(npas=1;npas<=maxpas;npas++)
	{
		//-- velocity verlet
		// modif_v_et_r(...)
		modif_v_et_r( pos, vit, acc);
		
		// eval_forces(...)
		V = eval_forces( force, acc, pos, &D );
		
		// modif_v_only(...)
		modif_v_only( vit, acc );
		K = eval_K( vit );
		
		T = C * K;
		
		//-- thermostat ?
		Real alpha=sqrt(K_th/K);

		if ( npas>=1500 && npas<1999)
		{			
			 T2 += T;
			 label2++;
		}
		else if ( npas==1999 ) //moyenne auto sur les T pour calculer la T de thermalisation 
		{
			T2  = T2 / label2;
			a   = (0.005-T2)/2500;  //coef de la droite de decroissance de T
			b   =  0.005 - 4500* a; // a*x + b
		}	
		else if( npas>=2000 && npas<4500 )
		{
			T1 = b + npas * a;			
			for(i=0;i<nat;i++)
			{
				vit[i].x = vit[i].x * alpha;
				vit[i].y = vit[i].y * alpha;
				vit[i].z = vit[i].z * alpha;
			}
			K    = K_th;
			T    = T1;       
			K_th = (3.0/2.0)*nat*kb*T1;
		}
		else if (npas>=4500)
		{
			K = K_th;
			T = T1;
		}
		
		P = (1.0 /(3.0*Vol)) * (K + D);
		
		//impression dans les fichier
		fichier_temp_press = fopen("temppress.txt","ab");   //pk dans temp.txt premier pas a 0??
		if(fichier_temp_press)
		{
			fprintf(fichier_temp_press,"%d  %lf  %lf\n",npas,T,P);    // fichier de temperature
			fclose(fichier_temp_press);
		}
		
		printf("npas= %d  E=%10g  V=%10g  K=%10g\n", npas, V+K, V, K);
		
		if( (npas%4) == 0 )
		{
			fichier = fopen("nrj.txt", "ab" );
			if( fichier )
			{
				fprintf( fichier, "%d %g %g %g\n",npas, V, K, V+K );
				fclose( fichier );
			}
			
			
			fichier = fopen("atoms.xyz", "ab");
			if(fichier)
			{
				fprintf( fichier, "%d\n", nat );
				fprintf( fichier, "\n" );
				for( i=0; i < nat; ++i )
				{
					const point *r = &pos[i];
					fprintf(fichier,"H %lf %lf %lf\n",r->x,r->y,r->z);
				}
				fclose(fichier);
			}
		}
		//break;
	}   //fin de la boucle principale
	
	
	
	free( force );
	free( acc );
	free( vit );
	free( pos );
	
	return 0;
}
